﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ExpressionPlotterControl;

namespace TubesMetnum4
{
    class Metode
    {
        public static double Trapesium(Expression func, double a, double b, int n)
        {
            double hasil = 0;
            double x = a;
            double selang = (b-a) / n;
            while (x < b)
            {
                if (!(Double.IsInfinity(func.Evaluate(x)) || Double.IsNaN(func.Evaluate(x))))

                hasil += selang * (func.Evaluate(x) + func.Evaluate(x + selang)) / 2;
                else
                    hasil += selang * (func.Evaluate(x + 0.000000001) + func.Evaluate(x + selang)) / 2;
                x += selang;
            }
            return hasil;
        }

        public static double Trapesium2(Func<double, double> func, double a, double b, int n)
        {
            double hasil = 0;
            double x = a;
            double selang = (b - a) / n;
            while (x < b)
            {
                if (!(Double.IsInfinity(func(x)) || Double.IsNaN(func(x))))
                    hasil += selang * (func(x) + func(x + selang)) / 2;
                else
                    hasil += selang * (func(x + 0.000000001) + func(x + selang)) / 2;
                x += selang;
            }
            return hasil;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="func">Fungsi yang ingin diselesaikan</param>
        /// <param name="a">nilai awal</param>
        /// <param name="b">nilai akhir</param>
        /// <param name="n">banyaknya interval yang diinginkan</param>
        /// <returns>Nilai integralnya</returns>
        public static double Simpson13(Expression func, double a, double b, int n)
        {
            double x = 0, F = 0;
            double h = (b - a) / n;
            double S = func.Evaluate(a);

            for (int i = 1; i < n; i += 2)
            {
                x = a + h * i;
                S += 4 * func.Evaluate(x);
            }

            for (int i = 2; i < n - 1; i += 2)
            {
                x = a + h * i;
                S += 2 * func.Evaluate(x);
            }

            S += func.Evaluate(b);
            F = h * S / 3;

            return F;
        }

        public static double Simpson132(Func<double, double> func, double a, double b, int n)
        {
            double x = 0, F = 0;
            double h = (b - a) / n;
            double S = func(a);

            for (int i = 1; i < n; i += 2)
            {
                x = a + h * i;
                S += 4 * func(x);
            }

            for (int i = 2; i < n - 1; i += 2)
            {
                x = a + h * i;
                S += 2 * func(x);
            }

            S += func(b);
            F = h * S / 3;

            return F;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="func">Fungsi yang ingin diselesaikan</param>
        /// <param name="a">nilai awal</param>
        /// <param name="b">nilai akhir</param>
        /// <param name="n">banyaknya interval yang diinginkan</param>
        /// <returns>Nilai integralnya</returns>
        public static double Simpson382(Func<double, double> func, double a, double b, int n)
        {
            double h, x, sigma, I;
            int r;

            h = (b - a) / n;
            x = a;
            I = func(a) + func(b);
            sigma = 0;
            for (r = 1; r < n - 1; r++)
            {
                x = x + h;
                if (r % 3 == 1)
                {
                    sigma = sigma + 2 * func(x);
                }
                else
                {
                    sigma = sigma + 3 * func(x);
                }
            }

            I = (I + sigma) * 3 * h / 8;
            return I;
        }
        public static double Simpson38(Expression func, double a, double b, int n)
        {
            double h, x, sigma, I;
            int r;

            h = (b - a) / n;
            x = a;
            I = func.Evaluate(a) + func.Evaluate(b);
            sigma = 0;
            for (r = 1; r < n - 1; r++)
            {
                x = x + h;
                if (r % 3 == 1)
                {
                    sigma = sigma + 2 * func.Evaluate(x);
                }
                else
                {
                    sigma = sigma + 3 * func.Evaluate(x);
                }
            }

            I = (I + sigma) * 3 * h / 8;
            return I;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="func">Fungsi yang ingin diselesaikan</param>
        /// <param name="a">nilai awal</param>
        /// <param name="b">nilai akhir</param>
        /// <param name="n">banyaknya interval yang diinginkan</param>
        /// <returns>Nilai integralnya</returns>
        public static double Roomberg(Expression func, double a, double b, int n)
        {
            double[] s = new double[n];
            int i,k;
            double var = 0;
            for (i = 1; i< n; i++)
                s[i] = 1;

            for (k = 1; k < n; k++)
            {
                for (i = 1; i <= k; i++)
                {
                    if (i == 1)
                    {
                        var = s[i];
                        s[i] = Trapesium(func, a, b, (int)Math.Pow(2, k - 1));     // sub-routine trapeze 
                    }                                       // integrated from 0 and 1
                    /* pow() is the number of subdivisions*/
                    else
                    {
                        s[k] = (Math.Pow(4, i - 1) * s[i - 1] - var) / (Math.Pow(4, i - 1) - 1);

                        var = s[i];
                        s[i] = s[k];
                    }
                }
            }
 
            return s[n - 1];
        }
        public static double Roomberg2(Func<double, double> func, double a, double b, int n)
        {
            double[] s = new double[n];
            int i, k;
            double var = 0;
            for (i = 1; i < n; i++)
                s[i] = 1;

            for (k = 1; k < n; k++)
            {
                for (i = 1; i <= k; i++)
                {
                    if (i == 1)
                    {
                        var = s[i];
                        s[i] = Trapesium2(func, a, b, (int)Math.Pow(2, k - 1));     // sub-routine trapeze 
                    }                                       // integrated from 0 and 1
                    /* pow() is the number of subdivisions*/
                    else
                    {
                        s[k] = (Math.Pow(4, i - 1) * s[i - 1] - var) / (Math.Pow(4, i - 1) - 1);

                        var = s[i];
                        s[i] = s[k];
                    }
                }
            }

            return s[n - 1];
        }

        /// <summary>
        /// Menghitung nilai integral suatu fungsi dengan kaidah titik tengah, galat O(h^2)
        /// </summary>
        /// <param name="func">Fungsi yang akan dihitung nilai integralnya</param>
        /// <param name="a">Batas awal</param>
        /// <param name="b">Batas akhir</param>
        /// <param name="n">Banyak pias</param>
        /// <returns></returns>
        public static double TitikTengah(Expression func, double a, double b, int n)
        {
            double h = (b - a) / n;
            double x = a + (h / 2);
            double sigma = func.Evaluate(x);

            for (int r = 1; r < n; r++)
            {
                x = x + h;
                sigma = sigma + func.Evaluate(x);
            }

            return (sigma * h);
        }
        public static double TitikTengah2(Func<double, double> func, double a, double b, int n)
        {
            double h = (b - a) / n;
            double x = a + (h / 2);
            double sigma = func(x);

            for (int r = 1; r < n; r++)
            {
                x = x + h;
                sigma = sigma + func(x);
            }

            return (sigma * h);
        }
        /// <summary>
        /// Menghitung nilai integrasi suatu fungsi dengan kaidah Gauss-Legendre
        /// </summary>
        /// <param name="func"></param>
        /// <param name="a">batas awal</param>
        /// <param name="b">batas akhir</param>
        /// <param name="n">orde metodeGauss-Legendre (2,3,4,atau 5)</param>
        /// <returns></returns>
        public static double GaussLegendre(Expression func, double a, double b, int n)
        {
            //transformasi selang
            Func<double, double> x = t => ((a + b) + (b - a) * t) / 2;
            double dx = (b - a) / 2;        //dx = (b-a)/2 dt
            double F = 0;
            
            switch (n)
            {
                case 2 :
                    F = dx * func.Evaluate(x.Invoke(1 / Math.Pow(3, 0.5))) + dx * func.Evaluate(x.Invoke(-1 / Math.Pow(3, 0.5)));
                    break;
                case 3 :
                    F = 5 / 9.0 * dx * func.Evaluate(x.Invoke(-1 * Math.Pow(0.6, 0.5))) + 8 / 9.0 * dx * func.Evaluate(x.Invoke(0)) + 5 / 9.0 * dx * func.Evaluate(x.Invoke(Math.Pow(0.6, 0.5)));
                    break;
                case 4:
                    F = 0.347854845 * dx * func.Evaluate(x.Invoke(-0.861136312)) + 0.652145155 * dx * func.Evaluate(x.Invoke(-0.339981044)) +
                         0.652145155 * dx * func.Evaluate(x.Invoke(0.339981044)) + 0.347854845 * dx * func.Evaluate(x.Invoke(0.861136312));
                    break;
                case 5:
                    F = 0.236926885 * dx * func.Evaluate(x.Invoke(-0.906179846)) + 0.478628670 * dx * func.Evaluate(x.Invoke(-0.538469310)) + 0.568888889 * dx * func.Evaluate(x.Invoke(0)) +
                         0.478628670 * dx * func.Evaluate(x.Invoke(0.538469310)) + 0.236926885 * dx * func.Evaluate(x.Invoke(0.906179846));
                    break;
                default :
                    F = double.MinValue;
                    break;
            }
            
            return F;
        }

        public static double GaussLegendre2(Func<double,double> func, double a, double b, int n)
        {
            //transformasi selang
            Func<double, double> x = t => ((a + b) + (b - a) * t) / 2;
            double dx = (b - a) / 2;        //dx = (b-a)/2 dt
            double F = 0;

            switch (n)
            {
                case 2:
                    F = dx * func(x.Invoke(1 / Math.Pow(3, 0.5))) + dx * func(x.Invoke(-1 / Math.Pow(3, 0.5)));
                    break;
                case 3:
                    F = 5 / 9.0 * dx * func.Invoke(x.Invoke(-1 * Math.Pow(0.6, 0.5))) + 8 / 9.0 * dx * func.Invoke(x.Invoke(0)) + 5 / 9.0 * dx * func.Invoke(x.Invoke(Math.Pow(0.6, 0.5)));
                    break;
                case 4:
                    F = 0.347854845 * dx * func(x.Invoke(-0.861136312)) + 0.652145155 * dx * func(x.Invoke(-0.339981044)) +
                         0.652145155 * dx * func(x.Invoke(0.339981044)) + 0.347854845 * dx * func(x.Invoke(0.861136312));
                    break;
                case 5:
                    F = 0.236926885 * dx * func(x.Invoke(-0.906179846)) + 0.478628670 * dx * func(x.Invoke(-0.538469310)) + 0.568888889 * dx * func(x.Invoke(0)) +
                         0.478628670 * dx * func(x.Invoke(0.538469310)) + 0.236926885 * dx * func(x.Invoke(0.906179846));
                    break;
                default:
                    F = double.MinValue;
                    break;
            }

            return F;
        }
    }
}
